# fig1-plot-garch.R
library(rugarch)
library(readstata13)
library(MASS) # for mvnorm
library(dynamac)
library(reshape2)
library(ggplot2)
library(gridExtra)
library(ggpubr)
setwd("~/Dropbox/Benton-Jordan-Philips/final-JOP-replication")


overall <- read.dta13("data/garch-estimates-Stata.dta")

names(overall)

estimated <- (overall$combination[!is.na(overall$estly)])


level.1 <- c("hetxatneg1", "hetxatzero", "hetxatpos1")
level.2 <- c("hetzatneg1", "hetzatzero", "hetzatpos1")
level.1.num <- level.2.num <- c(-1, 0, 1)

level.3 <- c("shockneg1", "shocknegpoint1", "shockpospoint1", "shockpos1")
level.3.num <- c(-1, -0.1, 0.1, 1)

the.paths <- rep(NA, length(level.1)*length(level.2)*length(level.3))*length(estimated)

counter <- 1
sim.data <- list()

for(i in 1:length(level.1)) {
    for(j in 1:length(level.2)) {
        for(k in 1:length(level.3)) {
            for(l in 1:length(estimated)) {
                the.paths[counter] <- paste0("data/figure1-data/", 
                                             paste(level.1[i], level.2[j], level.3[k], sep = "/"), "/comb_", estimated[l], ".dta")
                sim.data[counter][[1]] <- read.dta13(paste(the.paths[counter]))
                sim.data[counter][[1]]$garch <- na.omit(overall$garch[overall$combination == estimated[l]])
                sim.data[counter][[1]]$seed <- na.omit(overall$seed[overall$combination == estimated[l]])
                sim.data[counter][[1]]$combination <- na.omit(overall$combination[overall$combination == estimated[l]])
                sim.data[counter][[1]]$sd_of_hetx <- na.omit(overall$sd_of_hetx[overall$combination == estimated[l]])
                sim.data[counter][[1]]$sd_of_hetz <- na.omit(overall$sd_of_hetz[overall$combination == estimated[l]])
                sim.data[counter][[1]]$beta_hetx <- na.omit(overall$beta_hetx[overall$combination == estimated[l]])
                sim.data[counter][[1]]$beta_hetz <- na.omit(overall$beta_hetz[overall$combination == estimated[l]])
                sim.data[counter][[1]]$hetx_at <- level.1.num[i]
                sim.data[counter][[1]]$hetz_at <- level.2.num[j]
                sim.data[counter][[1]]$shock <- level.3.num[k]
                counter <- counter + 1
            }
        }
    }
}

all.sims.data <- data.frame(do.call("rbind", sim.data))
all.sims.data$GARCH <- all.sims.data$garch

all.sims.data.oneseed <- all.sims.data[all.sims.data$seed == unique(all.sims.data$seed)[2] & 
                                           all.sims.data$shock %in% c(-1, 1) & all.sims.data$GARCH %in% c(0.2, 0.5, 0.8) & 
                                           all.sims.data$sd_of_hetx %in% c(1) & all.sims.data$sd_of_hetz %in% c(1) &
                                           all.sims.data$beta_hetx %in% c(1) & all.sims.data$beta_hetz %in% c(1) &
                                           all.sims.data$time %in% seq(1, 10, 1),]

all.sims.data.oneseed$comb



all.sims.data.oneseed$ZAT <- factor(all.sims.data.oneseed$hetz_at, levels = c("-1", "0", "1"), labels = c("Z Held at Mean -1 SD", "Z Held at Mean", "Z Held at Mean +1 SD"))
all.sims.data.oneseed$XAT <- factor(all.sims.data.oneseed$hetx_at, levels = c("-1", "0", "1"), labels = c("X Held at Mean -1 SD", "X Held at Mean", "X Held at Mean +1 SD"))
all.sims.data.oneseed$GARCHlab <- factor(all.sims.data.oneseed$GARCH, levels = c("0.2", "0.5", "0.8"), labels = c("GARCH = 0.2", "GARCH = 0.5", "GARCH = 0.8"))
all.sims.data.oneseed$SHOCK <- factor(all.sims.data.oneseed$shock, levels = c("1", "-1"), labels = c("+1 SD", "-1 SD"))

#pdf("/Users/scj0014/Myfiles/Dropbox/Benton-Jordan-Philips/Picture Perfect/text/simulated.pdf", width = 7, height = 5)
ggplot(data = all.sims.data.oneseed, aes(x = time, y = sigma2_, color = as.factor(XAT), lty = as.factor(SHOCK))) + 
    geom_line(size = 1) + 
    scale_colour_grey(start = 0.7, end = 0) +
    facet_grid(vars(GARCHlab), vars(ZAT), scales = "free") +
    geom_vline(xintercept = 3) +
    scale_x_continuous(breaks = c(1, 3, 5, 7, 9), labels = c(1, 3, 5, 7, 9)) +
    theme(panel.grid.major.x = element_blank(), panel.grid.minor.x = element_blank(), 
          panel.background = element_rect(fill = "white", colour = NA), 
          panel.border = element_rect(fill = NA, colour = "grey20"), 
          panel.grid = element_line(colour = "grey92"), 
          panel.grid.minor = element_line(size = rel(0.5)), 
          strip.background = element_rect(fill = "grey85", colour = "grey20"), 
          legend.key = element_rect(fill = "white", colour = NA), complete = TRUE) +	
    labs(color = "X Value", lty = "Value of Shock to X", x = "Time Period 1 to 10 (Shock at t = 3)", y = "Conditional Variance")
dev.off()


